% matlab example of VAR
% > Econometric Toolbox > Model Selection > Specification Testing

% 1. load data
% clc; clear all;
% [val txt] = xlsread('../data/data_ls_macro.xlsx',2); save '../data/data2';
% fprintf('saved data');

%%
clc; clear all;load '../data/data2';
GDP = val(:,3); %real GDP
M2 = val(:,4);
Credit = val(:,5);
CPI = val(:,7);
lsLNH = val(:,2);

dates = datenum(txt(2:end,1));

Y = [Credit./CPI,M2./CPI, CPI, lsLNH];
X = dates;

% Transforming Data for Stationarity
% plot the data for trend
figure
subplot(2,2,1)
plot(dates,Y(:,1),'r');
title('Credit')
datetick('x')
grid on
subplot(2,2,2);
plot(dates,Y(:,2),'b');
title('M2')
datetick('x')
grid on
subplot(2,2,3);
plot(dates, Y(:,3), 'k')
title('CPI')
datetick('x')
grid on
subplot(2,2,4);
plot(dates, Y(:,4), 'magenta')
title('interbank interest rate 1M')
datetick('x')
grid on
hold off


% transfromed and plot again
Y = [diff(log(Y(:,1:2))), diff(Y(:,3)) diff(Y(1:end,4))];
X = dates(2:end);

figure
subplot(2,2,1)
plot(X,Y(:,1),'r');
title('Credit')
datetick('x')
grid on
subplot(2,2,2);
plot(X,Y(:,2),'b');
title('M2')
datetick('x')
grid on
subplot(2,2,3);
plot(X, Y(:,3),'k'),
title('CPI')
datetick('x')
grid on
subplot(2,2,4);
plot(X, Y(:,4),'magenta'),
title('interbank 1M')
datetick('x')
grid on


% plot in one figure only

Y(:,1:2) = 100*Y(:,1:2);
figure
plot(X,Y(:,1),'r');
hold on
plot(X,Y(:,2),'b');
datetick('x')
grid on
plot(X,Y(:,3),'k');
plot(X,Y(:,4),'magenta');
legend('Credit','M2','CPI','interbank 1M');
hold off

%% check stationary
temp = diff(log(lsLNH));
% temp = diff(log(M2./CPI));
% temp = diff(CPI);
fdates = size(dates,1) - size(temp,1);
x = dates(fdates+1:end);


close all;plot(x,temp);datetick('x',27);axis tight;

[x1, pvalue] = adftest(temp); [x1, pvalue]
% -> result :[ lsLNH CPI, log(Credit/CPI), log(M2/CPI) ] -> I(1)
%           : GDP (adjusted and stationary) -> I(2)


%% selecting and fitting data

% === SELECTING MODEL
% make the series same length
dCredit = 100*diff(log(Credit./CPI));
dM2 = 100*diff(log(M2./CPI));
dCPI = diff(CPI);
dlsLNH = diff(lsLNH);
Y = [dCredit dM2 dCPI dlsLNH];



%create 4 model
dt = logical(eye(4));
VAR2diag = vgxset('ARsolve',repmat({dt},2,1),...
    'asolve',true(4,1),'Series',{'Credit','M2','CPI','1-mo interbank'},'nAR',2);
VAR2full = vgxset(VAR2diag,'ARsolve',[]);
VAR4diag = vgxset(VAR2diag,'nAR',4,'ARsolve',repmat({dt},4,1));
VAR4full = vgxset(VAR2full,'nAR',4);
VAR3full = vgxset(VAR2full,'nAR',3);

% === Choosing Presample, Estimation, and Forecast Periods. 

Ypre = Y(1:4,:);
T = floor(.9*size(Y,1));
Yest = Y(5:T,:);
YF = Y((T+1):end,:);
TF = size(YF,1);

% fitting
[EstSpec1,EstStdErrors1,LLF1,W1] = ...
    vgxvarx(VAR2diag,Yest,[],Ypre,'CovarType','Diagonal');
[EstSpec2,EstStdErrors2,LLF2,W2] = ...
    vgxvarx(VAR2full,Yest,[],Ypre);
[EstSpec3,EstStdErrors3,LLF3,W3] = ...
    vgxvarx(VAR4diag,Yest,[],Ypre,'CovarType','Diagonal');
[EstSpec4,EstStdErrors4,LLF4,W4] = ...
    vgxvarx(VAR4full,Yest,[],Ypre);
[EstSpec5,EstStdErrors5,LLF5,W5] = ...
    vgxvarx(VAR3full,Yest,[],Ypre);

%% Checking Model Adequacy

[isStable1,isInvertible1] = vgxqual(EstSpec1);
[isStable2,isInvertible2] = vgxqual(EstSpec2);
[isStable3,isInvertible3] = vgxqual(EstSpec3);
[isStable4,isInvertible4] = vgxqual(EstSpec4);
[isStable5,isInvertible5] = vgxqual(EstSpec5);
[isStable1,isStable2,isStable3,isStable4, isStable5]

% Likelihood Ratio Tests
% count the active number parameter
[n1,n1p] = vgxcount(EstSpec1);
[n2,n2p] = vgxcount(EstSpec2);
[n3,n3p] = vgxcount(EstSpec3);
[n4,n4p] = vgxcount(EstSpec4);
[n5,n5p] = vgxcount(EstSpec5);
% test
reject1 = lratiotest(LLF2,LLF1,n2p - n1p)
reject3 = lratiotest(LLF4,LLF3,n4p - n3p)
reject4 = lratiotest(LLF4,LLF2,n4p - n2p)
reject5 = lratiotest(LLF5,LLF2,n5p - n2p)


% AIC 
AIC = aicbic([LLF1 LLF2 LLF3 LLF4 LLF5],[n1p n2p n3p n4p n5p])

%%
% Compare forecast with forecast period data
Yest = Y(5:T,:);
[FY1,FYCov1] = vgxpred(EstSpec1,TF,[],Yest);
[FY2,FYCov2] = vgxpred(EstSpec2,TF,[],Yest);
[FY3,FYCov3] = vgxpred(EstSpec3,TF,[],Yest);
[FY4,FYCov4] = vgxpred(EstSpec4,TF,[],Yest);
[FY5,FYCov5] = vgxpred(EstSpec5,TF,[],Yest);

%% rolling forecasting technique
numPredictors = size(Y,2);
FY1 = zeros(TF,numPredictors);
FY2 = zeros(TF,numPredictors);
FY3 = zeros(TF,numPredictors);
FY4 = zeros(TF,numPredictors);
FY5 = zeros(TF,numPredictors);
Yest = Y(5:T,:);
for i = 1:TF
    FY1(i,:) = vgxpred(EstSpec1,1,[],Yest);
    FY2(i,:) = vgxpred(EstSpec2,1,[],Yest);
    FY3(i,:) = vgxpred(EstSpec3,1,[],Yest);
    FY4(i,:) = vgxpred(EstSpec4,1,[],Yest);
    FY5(i,:) = vgxpred(EstSpec5,1,[],Yest);
    
    if(i~=TF)
        Yest = Y(5:T+i,:);
        %estimate again
        [EstSpec1,EstStdErrors1,LLF1,W1] = ...
            vgxvarx(VAR2diag,Yest,[],Ypre,'CovarType','Diagonal');
        [EstSpec2,EstStdErrors2,LLF2,W2] = ...
            vgxvarx(VAR2full,Yest,[],Ypre);
        [EstSpec3,EstStdErrors3,LLF3,W3] = ...
            vgxvarx(VAR4diag,Yest,[],Ypre,'CovarType','Diagonal');
        [EstSpec4,EstStdErrors4,LLF4,W4] = ...
            vgxvarx(VAR4full,Yest,[],Ypre);
        [EstSpec5,EstStdErrors5,LLF5,W5] = ...
            vgxvarx(VAR3full,Yest,[],Ypre);
    end
end

%%

figure
% vgxplot(EstSpec2,Yest,FY2,FYCov2)
% vgxplot(EstSpec4,Yest,FY4,FYCov4)
% Sum Square Error
error1 = YF - FY1;
error2 = YF - FY2;
error3 = YF - FY3;
error4 = YF - FY4;
error5 = YF - FY5;

SSerror1 = error1(:)' * error1(:);
SSerror2 = error2(:)' * error2(:);
SSerror3 = error3(:)' * error3(:);
SSerror4 = error4(:)' * error4(:);
SSerror5 = error5(:)' * error5(:);
figure
% bar([SSerror1 SSerror2 SSerror3 SSerror4 SSerror5],.5)
bar([SSerror2 SSerror4 SSerror5],.5)
ylabel('Sum of squared errors')
set(gca,'XTickLabel',...
    {'AR2 full' 'AR4 full','AR3 full'})
title('Sum of Squared Forecast Errors')

vgxdisp(EstSpec2)

%% % choosing lag between 1-12:
% Ybackup = Y;
% Y = Y(:,[1,2,4]);

T = floor(.9*size(Y,1));

Ypre1 = Y(1:12,:);
Yest1 = Y(13:T,:);
result = zeros(12,4); %likelihood and number active parameters
YF = Y((T+1):end,:);
TF = size(YF,1);
numPredictors = size(Y,2);

for iLag = 1: 12
    VARtemp = vgxset('asolve',true(numPredictors,1),'nAR',iLag);
%     VARtemp = vgxset(VARtemp,'Series',{'Credit','M2','1-mo interbank'});
    [EstSpecTemp,EstStdErrorsTemp,LLFTemp,WTemp] = vgxvarx(VARtemp,Yest1,[],Ypre1);
    [isStable,isInvertible] = vgxqual(EstSpecTemp);
    [nTemp,nTempp] = vgxcount(EstSpecTemp);
    result(iLag,1:3) = [LLFTemp, nTempp,isStable];
    
    %calculate rolling forecast performance
    

    FY1 = zeros(TF,numPredictors);
    
    % numPredictors = size(Y,2);
    EstSpecTemp1 = EstSpecTemp;
    Yest2 = Y(13:T,:);
    for i = 1:TF
        FY1(i,:) = vgxpred(EstSpecTemp1,1,[],Yest2);        
        if(i~=TF)
            Yest2 = Y(13:T+i,:);
            %estimate again
            [EstSpecTemp1,EstStdErrorsTemp1,LLFTemp1,WTemp1] = ...
                vgxvarx(VARtemp,Yest2,[],Ypre1);
        end
    end
    error1 = YF - FY1;
    result(iLag,4) = error1(:)' * error1(:);
end
[AIC,BIC] = aicbic(result(:,1),result(:,2),size(Yest1,1));
HQC = 2.*log(log(size(Yest1,1))).*result(:,2) - 2.*result(:,1);
xxx = [AIC BIC HQC result(:,4)];

%%
close all;figure
plot(1:12,xxx(:,1:3)); hold on;
bar(1:12, xxx(:,4))
xlabel('Lag')
title('Information Criterion for VAR''s lag selection')
legend('AIC','BIC','HQC','error');

%% Forecasting
bestLag = 1;
finalModelSpec = vgxset('asolve',true(numPredictors,1),'nAR',bestLag);
finalModelSpec = vgxset(finalModelSpec,'Series',{'Credit','M2','CPI','1-mo interbank'});
[EstSpecFinal,EstStdErrorsFinal,LLFFinal,WFinal] = vgxvarx(finalModelSpec,Y(2:end,:),[],Y(1:1,:));
vgxdisp(EstSpecFinal);



%% forecast with no error bars
[ypred, ycov] = vgxpred(EstSpecFinal,10,[],YF);
y_raw = [Credit./CPI, M2./CPI, CPI, lsLNH];
yfirst = [log(y_raw(end,1:2)) ,y_raw(end,3),y_raw(end,4)];
ypred1 = [yfirst;ypred(:,1:2)/100 ypred(:,3:4)];
ypred2 = cumsum(ypred1);
ypred2(:,1:2) = exp(ypred2(:,1:2));
lasttime = dates(end);
timess = lasttime:31:lasttime+310;

figure;
for i = 1:4
    subplot(2,2,i);
    plot(dates,y_raw(:,i)); hold on
    plot(timess,ypred2(:,i),'r-');hold off
    title(finalModelSpec.Series{i});    
    legend('Data','Forecast','Location','NorthWest');
    %draw vertical line
    hx = graph2d.constantline(lasttime, 'LineStyle',':', 'Color',[.7 .7 .7]);
    changedependvar(hx,'x');
    datetick('x','mm-yy');
end


%% forecasting with error bars
rng(1);
ysim0 = vgxsim(EstSpecFinal,10,[],YF,[],5000);

y_raw = [Credit./CPI, M2./CPI, CPI, lsLNH];
yfirst = [log(y_raw(end,1:2)) ,y_raw(end,3),y_raw(end,4)];
ysim0(:,1:2,:) = ysim0(:,1:2,:)./100;

ysim = [repmat(yfirst,[1,1,5000]);ysim0];
ysim2 = cumsum(ysim(:,1:end,:));
ysim2(:,1:2,:) = exp(ysim2(:,1:2,:));

ymean = mean(ysim2,3);
ystd = std(ysim2,0,3);


figure;
for i = 1:4
    subplot(2,2,i);
    plot(dates,y_raw(:,i)); hold on
    plot(timess,ymean(:,i),'r-');
    plot(timess,ymean(:,i)-ystd(:,i),'g--');
    plot(timess,ymean(:,i)+ystd(:,i),'g--');
    hold off
    title(finalModelSpec.Series{i});    
%     legend('Data','Forecast','Location','NorthWest');
    %draw vertical line
    hx = graph2d.constantline(lasttime, 'LineStyle',':', 'Color',[.7 .7 .7]);
    changedependvar(hx,'x');
    datetick('x','mm-yy');
end

%% Calculating Impulse Responses
